library(foreign)
setwd("~/Google Drive/Broockman/Completed/Extremity vs Consistency/ReplicationData/PanelSurveyApril2014")

data <- read.dta("ssi-spatial-all-merged-with-wave1.dta")
names(data)

save <- TRUE

wave1 <- data[,12:23] - 4
wave2 <- data[,25:36] - 4
least.extreme.response <- matrix(ncol = ncol(wave1), nrow = nrow(wave1))
most.extreme.response <- matrix(ncol = ncol(wave1), nrow = nrow(wave1))
for(row in 1:nrow(wave1)){
	for(col in 1:ncol(wave1)){
		if(!is.na(wave1[row,col]) & !is.na(wave2[row,col])){
			if(abs(wave1[row,col]) >= abs(wave2[row,col])){
				least.extreme.response[row,col] <- wave2[row,col]
				most.extreme.response[row,col] <- wave1[row,col]	
			}
			else{
				least.extreme.response[row,col] <- wave1[row,col]
				most.extreme.response[row,col] <- wave2[row,col]
			}
		}
	}
}
wave1 <- wave1 + 4
wave2 <- wave2 + 4
least.extreme.response <- least.extreme.response + 4
most.extreme.response <- most.extreme.response + 4
mean.response <- (wave1 + wave2)/2

weight <- data$weight


#Regressions
pktest <- function(x){
	x <- abs(x-4)
	row.averages <- rowMeans(x, na.rm = TRUE)
	summary(lm(row.averages ~ data$pkscale, weights=data$weight))
}
pktest(wave1)
pktest(wave2)
pktest(least.extreme.response)
pktest(mean.response)

getColProps <- function(wave.col, array = c(1:7)){
	result <- vector("numeric", length=length(array))
	for(i in 1:length(array)){
		this.row <- as.numeric(wave.col==array[i])
		result[i] <- weighted.mean(this.row, weight, na.rm=TRUE)
	}
	return(result)
}
wave1.issue.freqs <- apply(wave1, 2, getColProps)
wave2.issue.freqs <- apply(wave2, 2, getColProps)
mebound.issue.freqs <- apply(least.extreme.response, 2, getColProps)
most.extreme.issue.freqs <- apply(most.extreme.response, 2, getColProps)



#Naming
issue.names <- c("Environment","Social Security","Guns","Health","Immigration","Taxes","Abortion","Medicare","Gay Rights","Unions","Contraception","Education")

#Mean response across two waves
array <- seq(from=1,to=7,by=.5)
mean.response.issue.freqs <- apply(mean.response, 2, getColProps, array=array)

#Figures
makeplot13 <- function(thirteen.point, title, ylimtop=.3, sub="", cex=.6){
	bp <- barplot(thirteen.point, main = title, ylim=c(0,ylimtop), names.arg = array, col = c("blue4", "blue3", "blue2", "blue1", "blue", "blue", "Purple", "red", "red", "red1", "red2", "red3", "red4"), sub=sub)
	text(bp, thirteen.point + .01, paste0(round(thirteen.point*100, 0),"%"),cex=cex,pos=3)
}

if(save) pdf("average_response_weighted.pdf", width = 900*2/50, height = 580*2/50, pointsize = 16*2)
par(mfrow=c(3,4), oma = c(0,0,3,0), mar = c(3,3,3,3))
for(i in 1:ncol(mean.response.issue.freqs)){
	makeplot13(mean.response.issue.freqs[,i], issue.names[i])
}
par(mfrow=c(1,1), oma = c(0,0,0,0))
title(main="Average of Wave 1 and Wave 2 Responses")
if(save) dev.off()


if(save) pdf("average_response_average_issue.pdf", width = 700*2/50, height = 500*2/50, pointsize = 16*2)
makeplot13(rowMeans(mean.response.issue.freqs), "Distribution of Responses on 'Average Issue' - Mass Public", sub="Mean of Wave 1 and Wave 2 Responses", ylimtop = .2, cex=.9)
if(save) dev.off()


makeplot <- function(seven.point, title, ylimtop = .45, sub=""){
	bp <- barplot(seven.point, main = title, ylim=c(0,ylimtop), names.arg = c(1:7), col = c("darkblue", "blue3", "blue", "Purple", "red", "red3", "darkred"), sub=sub)
	text(bp, seven.point + .01, paste0(round(seven.point*100, 0),"%"),cex=.9,pos=3)
}


if(save) pdf("dist_of_multis_weighted.pdf", width = 900*2/50, height = 580*2/50, pointsize = 16*2)
par(mfrow=c(3,4), oma = c(0,0,3,0), mar = c(3,3,3,3))
for(i in 1:ncol(wave1.issue.freqs)){
	makeplot(wave1.issue.freqs[,i], issue.names[i])
}
par(mfrow=c(1,1), oma = c(0,0,0,0))
title(main="Distribution of Opinion on Twelve Issues - Mass Public")
if(save) dev.off()

if(save) pdf("mass_mebound_issues_then_respondents.pdf", width = 700*2/50, height = 500*2/50, pointsize = 16*2)
makeplot(rowMeans(mebound.issue.freqs), "Distribution of Responses on 'Average Issue' - Mass Public", sub="Measurement Error Model Lower Bound", ylimtop = .4)
if(save) dev.off()

if(save) pdf("most_extreme_response_issues_then_respondents.pdf", width = 700*2/50, height = 500*2/50, pointsize = 16*2)
makeplot(rowMeans(most.extreme.issue.freqs), "Distribution of Responses on 'Average Issue' - Mass Public", sub="Measurement Error Model Upper Bound", ylimtop = .4)
if(save) dev.off()

if(save) pdf("most_extreme_dist_of_multis_weighted.pdf", width = 900*2/50, height = 580*2/50, pointsize = 16*2)
par(mfrow=c(3,4), oma = c(0,0,3,0), mar = c(3,3,3,3))
for(i in 1:ncol(wave1.issue.freqs)){
	makeplot(most.extreme.issue.freqs[,i], issue.names[i], ylimtop = .5)
}
par(mfrow=c(1,1), oma = c(0,0,0,0))
title(main="Measurement Error Upper Bound - Mass Public")
if(save) dev.off()